library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.6 v dplyr 1.0.8
## v tidyr 1.2.0 v stringr 1.4.0
## v readr 2.1.2 v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(lubridate)
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(mapview)
library(maptools)
## Loading required package: sp
## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired by the end of 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
library(ggrepel)
library(patchwork)
library(rayshader)
library(osmdata)
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(ggmap)
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
Primero descargamos datos de los viajes realizados en el sistema de transporte metropolitano (STM) en Montevideo, como el archivo es grande lo descargamos y descomprimimos directamente de la pagina del catalogo de datos abiertos
temp <- tempfile()
download.file("https://catalogodatos.gub.uy/dataset/b1b22d81-9333-4a1b-8254-589268a698bf/resource/aa84ab90-7934-49ae-a330-57403f7e4e2e/download/viajes_stm_072022.zip",temp)
datos_julio <- read.csv(unz(temp, "viajes_stm_072022.csv"))
unlink(temp)
Cargamos Shapefiles con recorridos de los omnibus, paradas y Zonas censales para contextualizar geograficamente los datos previos dado que no se cuenta en los mismos con un x e y
zonas_censales <- st_read("data/censo/ine_seg_11.shp")
## Reading layer `ine_seg_11' from data source
## `E:\Gonzalo Machin\Maestria\Materias\CDC 2\tps_cdc_2\data\censo\ine_seg_11.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 4313 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 366582.2 ymin: 6127919 xmax: 858252.1 ymax: 6671738
## CRS: NA
paradas <- st_read("data/paradas_shp/v_uptu_paradas.shp", options = "ENCODING=UTF-8")
## options: ENCODING=UTF-8
## Reading layer `v_uptu_paradas' from data source
## `E:\Gonzalo Machin\Maestria\Materias\CDC 2\tps_cdc_2\data\paradas_shp\v_uptu_paradas.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 38764 features and 10 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 554395.5 ymin: 6134611 xmax: 591653.9 ymax: 6158038
## Projected CRS: WGS 84 / UTM zone 21S
recorridos <- st_read("data/recorridos/v_uptu_lsv.shp")
## Reading layer `v_uptu_lsv' from data source
## `E:\Gonzalo Machin\Maestria\Materias\CDC 2\tps_cdc_2\data\recorridos\v_uptu_lsv.shp'
## using driver `ESRI Shapefile'
## replacing null geometries with empty geometries
## Simple feature collection with 652 features and 8 fields (with 2 geometries empty)
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: 554322.5 ymin: 6134518 xmax: 591684.6 ymax: 6158125
## Projected CRS: WGS 84 / UTM zone 21S
recorridos_no_maximales <- st_read("data/recorridos_no_maximales/uptu_variante_no_maximal.shp")
## Reading layer `uptu_variante_no_maximal' from data source
## `E:\Gonzalo Machin\Maestria\Materias\CDC 2\tps_cdc_2\data\recorridos_no_maximales\uptu_variante_no_maximal.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 710 features and 13 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: 555656.7 ymin: 6134558 xmax: 591668.1 ymax: 6158125
## Projected CRS: WGS 84 / UTM zone 21S
Corregimos crs y cambiamos todos para trabajar unicamente en CRS = 4326
zonas_censales <- zonas_censales %>%
st_set_crs(32721) %>%
st_transform(4326)
paradas <- paradas %>%
st_transform(4326)
recorridos <- recorridos %>%
st_transform(4326)
recorridos_no_maximales <- recorridos_no_maximales %>%
st_transform(4326)
Filtramos los datos de las zonas censales para unicamente trabajar con los datos del departamento de Montevideo
zonas_censales <- zonas_censales %>%
filter(NOMBDEPTO == "MONTEVIDEO")
Observamos los datos geograficos al momento
ggplot()+
geom_sf(data = zonas_censales)+
geom_sf(data = recorridos)+
geom_sf(data = recorridos_no_maximales)+
geom_sf(data = paradas)+
theme_bw()
Observamos que el STM tiene alcance mas alla de las fronteras del departamento, pero en este trabajo unicamente utilizaremos aquellas paradas que estan dentro del departamento, para esto hacemos un spatial join para poder filtrar las paradas que no estan en el departamento y tambien asignarles el codigo de la zona censal dentro de la que pertenecen, asi podemos empezar a trabajar con los datos de transporte publico dentro de cada zona censal.
sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
paradas <- paradas %>%
st_make_valid() %>%
st_join(zonas_censales)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
paradas <- paradas %>%
filter(!is.na(CODSEG))
Observamos los datos de las secciones censales para ver en que nivel de detalle analizamos su subdivision
head(zonas_censales)
## Simple feature collection with 6 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -56.20501 ymin: -34.91282 xmax: -56.12181 ymax: -34.89936
## Geodetic CRS: WGS 84
## AREA PERIMETER DEPTO SECCION SEGMENTO LOCALIDAD CODSEC CODSEG CODLOC
## 1 14476.50 1785.1361 1 0 0 20 100 100000 1020
## 2 122200.01 1563.4714 1 1 1 20 101 101001 1020
## 3 151608.67 1676.3479 1 1 2 20 101 101002 1020
## 4 99377.28 1456.9389 1 1 3 20 101 101003 1020
## 5 54395.51 1040.6595 1 1 104 20 101 101104 1020
## 6 41055.69 838.3599 1 1 105 20 101 101105 1020
## NOMBDEPTO NOMBLOC CDEPTO_ISO CLOC_ISO geometry
## 1 MONTEVIDEO MONTEVIDEO UYMO UYMOMON MULTIPOLYGON (((-56.13303 -...
## 2 MONTEVIDEO MONTEVIDEO UYMO UYMOMON MULTIPOLYGON (((-56.20449 -...
## 3 MONTEVIDEO MONTEVIDEO UYMO UYMOMON MULTIPOLYGON (((-56.1998 -3...
## 4 MONTEVIDEO MONTEVIDEO UYMO UYMOMON MULTIPOLYGON (((-56.19956 -...
## 5 MONTEVIDEO MONTEVIDEO UYMO UYMOMON MULTIPOLYGON (((-56.19834 -...
## 6 MONTEVIDEO MONTEVIDEO UYMO UYMOMON MULTIPOLYGON (((-56.19813 -...
zonas_censales <- zonas_censales %>%
mutate(CODSEC = as.factor(CODSEC))
ggplot()+
geom_sf(data = zonas_censales,aes(fill=CODSEC))+
labs(title = "Secciones Departamentles",
y="Latitud",
x="Longitud",
legend="Seccion",
caption="Fuente: INE")+
theme_bw()
Observamos que tiene 3 niveles de subdivision, el departamento, la seccion y el segmento. como el departamento solo usaremos uno observaremos el nivel de detalle de cada seccion para observar si con este nivel es suficiente para sacar conclusiones geograficas o es necesario ir un nivel mas abajo y analiza por segmento censal
zonas_censales_CODSEC <- zonas_censales %>%
group_by(CODSEC) %>%
summarise()
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
ggplot()+
geom_sf(data = zonas_censales_CODSEC,aes(fill=CODSEC))+
geom_rect(aes(xmin = -56.235, xmax = -56.11, ymin = -34.85, ymax = -34.94), color = "red", fill = NA)+
labs(title = "Secciones Departamentles",
y="Latitud",
x="Longitud",
legend="Seccion",
caption="Fuente: INE")+
theme_bw()
Para facilitar el calculo y el manejo de 22 millones de datos utilizaremos group by y summarise para resumir el dataset y poder graficar los datos por zona censal
datos_julio_resumido <- datos_julio %>%
group_by(codigo_parada_origen,
sevar_codigo,
dsc_linea) %>%
summarise(Total = sum(cantidad_pasajeros)) %>%
rename(COD_UBIC_P = codigo_parada_origen) %>%
mutate(parada_linea = paste(dsc_linea,sevar_codigo,COD_UBIC_P,sep="_")) %>%
select(parada_linea,
Total)
## `summarise()` has grouped output by 'codigo_parada_origen', 'sevar_codigo'. You
## can override using the `.groups` argument.
## Adding missing grouping variables: `COD_UBIC_P`, `sevar_codigo`
paradas <- paradas %>%
mutate(parada_linea = paste(DESC_LINEA,COD_VARIAN,COD_UBIC_P,sep="_"))
paradas <- paradas %>%
left_join(datos_julio_resumido,by = "parada_linea") %>%
mutate_each(funs(replace(., which(is.na(.)), 0))) %>%
mutate(CODSEC = as.character(CODSEC))
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
## Warning: `mutate_each_()` was deprecated in dplyr 0.7.0.
## Please use `across()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
ggplot(paradas %>%
st_set_geometry(NULL) %>%
group_by(CODSEC) %>%
summarise(Total = sum(Total)))+
geom_bar(aes(y=reorder(CODSEC,(Total)),weight = Total),fill = "#60E6C8")+
labs(title="Total de Viajes Vendidos por Seccion Departamental",
y="Seccion",
x="Total",
caption="Fuente: STM")+
theme_bw()
Observamos que las secciones 110 y 111 son las dos que mas viajes vendidos tienen, ahora observamos espacialmente
viajes_por_seccion <- paradas %>%
st_set_geometry(NULL) %>%
group_by(CODSEC) %>%
summarise(Total = sum(Total))
zonas_censales_CODSEC <- left_join(zonas_censales_CODSEC,viajes_por_seccion,by = "CODSEC")
ggplot()+
geom_sf(data = zonas_censales_CODSEC,aes(fill=Total))+
scale_fill_gradient(low = "#FDFF62",high = "#FF6262")+
labs(title = "Total de Viajes Vendidos por Zona Censal",
x = "Longitud",
y = "Latitud",
legend = "Viajes Vendidos",
caption = "Fuente: STM")+
theme_bw()
Se puede observar que la zona sureste de la ciudad es la seccion con mas viajes vendidos, superando los dos millones de viajes en el mes de Julio 2022, para concluir aun mas sobre estas zonas se puede hacer mas hincapie en relacion a la poblacion y puestos de trabajo en cada zona censal, pero sobre esto se indagara en los proximos trabajos practicos.
En la segunda entraga de este trabajo practico abordaremos que tan bueno es el alcance del transporte publico y entenderemos si todos los rincones del departamento tienen la misma accesibilidad al sistema
Para esto calcularemos la distancia a la parada mas cercana de diferentes puntos de las zonas censales
Primero responderemos una pregunta que quedo pendiente del TP anterior, existe una relacion entre la cantidad de viajes vendidos y la densidad de poblacion o la superficie?
Para esto necesitaremos los datos del censo 2011
Poblacion <- read.csv("data/censo/personas_por_zona.csv")
Con los datos de poblacion cargados, procedemos a limpiarlos un poco y hacer un join para asignarlos a cada seccion censal
Poblacion <- Poblacion %>%
mutate(CODSEG = ((DPTO*100000) + (SECC*1000) + SEGM)) %>%
mutate(CODSEC = ((DPTO*100) + SECC))
Poblacion_por_seccion <- Poblacion %>%
group_by(CODSEC) %>%
summarise(Total = sum(Total))
zonas_censales_CODSEC <- zonas_censales_CODSEC %>%
mutate(CODSEC = as.numeric(CODSEC)) %>%
left_join(Poblacion_por_seccion,by = "CODSEC")
zonas_censales_CODSEC <- zonas_censales_CODSEC %>%
rename(Viajes_totales = Total.x,
Poblacion = Total.y)
zonas_censales_CODSEC <- zonas_censales_CODSEC %>%
filter(!is.na(Poblacion))
Calculamos el area de cada zona censal y la densidad de poblacion
zonas_censales_CODSEC <- zonas_censales_CODSEC %>%
mutate(Area_km2 = units::set_units(st_area(zonas_censales_CODSEC), km^2)) %>%
mutate(Area_km2 = sub(" [km^2].*","",Area_km2)) %>%
mutate(Area_km2 = as.numeric(Area_km2)) %>%
mutate(Densidad = Poblacion/Area_km2)
pob <- ggplot(zonas_censales_CODSEC,aes(x=Poblacion,y=Viajes_totales))+
geom_point()+
geom_smooth(method=lm, se=FALSE, col='red', size=0.5)+
labs(title = "Viajes Totales, Poblacion y Densidad",
x="Poblacion",
y = "Viajes Vendidos")+
theme_bw()
density <- ggplot(zonas_censales_CODSEC,aes(x=Densidad,y=Viajes_totales))+
geom_point()+
geom_smooth(method=lm, se=FALSE, col='red', size=0.5)+
labs(x="Densidad de Poblacion",
y = "Viajes Vendidos")+
theme_bw()
pob + density
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
Se puede observar una mayor correlacion y relacion lineal de la poblacion con la cantidad de viajes vendidos que la densidad.
Ahora estudiaremos la accesibilidad del Transporte Publico en la ciudad, para esto cargaremos primero las direcciones oficiales de la ciudad
direcciones <- st_read("data/direcciones/v_mdg_accesos.shp", options = "ENCODING=UTF-8")
## options: ENCODING=UTF-8
## Reading layer `v_mdg_accesos' from data source
## `E:\Gonzalo Machin\Maestria\Materias\CDC 2\tps_cdc_2\data\direcciones\v_mdg_accesos.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 372443 features and 7 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 552746.8 ymin: 6134530 xmax: 588662.6 ymax: 6159661
## Projected CRS: WGS 84 / UTM zone 21S
Observamos el dataset cargado
ggplot()+
geom_sf(data = zonas_censales_CODSEC)+
geom_sf(data = direcciones)+
theme_bw()
Ahora volveremos a transformar a CRS 32721 el .shp de paradas
paradas <- paradas %>%
st_transform(32721)
Para estudiar la accesibilidad crearemos un buffer de 500 metros de cada parada para saber que tantas lineas y variantes tiene acceso una direccion particular, en este caso cabe destacar que nuestro dataset de paradas crea un punto de iguales coordenadas para cada variante de linea de omnibus que para en esa parada, por lo tanto estariamos calculando el total de variantes a las que se tiene acceso desde una direccion particular..
paradas_buffer <- st_buffer(paradas, 500)
ggplot()+
geom_sf(data = zonas_censales_CODSEC)+
geom_sf(data = paradas_buffer)+
geom_sf(data = paradas)+
theme_bw()
Con el buffer podemos saber cuantas direcciones son abastecidas y a cuantas lineas tienen acceso, dado que queremos conservar aquellas direcciones que no tienen acceso al STM a menos de 500 mts, utilizaremos st_join para unir ambos datasets
Primero uniremos los buffers para mejorar la velocidad de calculo del st_join
paradas_buffer_lineas <- paradas_buffer %>%
group_by(DESC_LINEA) %>%
summarise(Paradas = n_distinct(COD_UBIC_P.x))
ggplot()+
geom_sf(data = zonas_censales_CODSEC)+
geom_sf(data = paradas_buffer_lineas)+
geom_sf(data = paradas)+
theme_bw()
Accesibilidad <- st_join(direcciones,paradas_buffer_lineas)
Ahora calculamos el total de Lineas a las que una direccion particular puede acceder a menos de 500 mts
Accesibilidad <- Accesibilidad %>%
mutate(Direccion = paste(NOM_CALLE,NUM_PUERTA,sep=" ")) %>%
mutate_all(~replace(., is.na(.), 0)) %>%
mutate(Accesos = if_else(Paradas == 0,0,1))
Direcciones_accesibles <- Accesibilidad %>%
st_set_geometry(NULL) %>%
group_by(Direccion) %>%
summarise(Lineas_accesibles = sum(Accesos))
Ahora procedemos a unir y graficar el total de paradas y lineas accesibles a menos de 500 metros desde cada direccion de la ciudad
direcciones <- direcciones %>%
mutate(Direccion = paste(NOM_CALLE,NUM_PUERTA,sep=" ")) %>%
left_join(Direcciones_accesibles,by = "Direccion")
my_breaks = c(0,2,5,15,50,200,500)
ggplot()+
geom_sf(data = zonas_censales_CODSEC)+
geom_sf(data = paradas_buffer_lineas)+
geom_sf(data = direcciones,aes(color = Lineas_accesibles))+
labs(title = "Lineas y Paradas Accesibles desde cada Direccion",
x = "Longitud",
y = "Latitud",
legend = "Accesibilidad")+
scale_color_gradient(trans = "log",breaks = my_breaks,labels = my_breaks)+
theme_bw()
## Warning: Transformation introduced infinite values in discrete y-axis
Se puede apreciar que aquellas direcciones que se ubican cerca o sobre de las principales arterias de la ciudad tienen una mayor accesibilidad al STM, mientras que en gris se pueden observar aquellos lugares que estan por fuera, predominantemente ubicadas en el norte y oeste de la ciudad, esto se puede explicar probablemente por la poblacion y baja densidad de poblacion.
Para observar mejor el impacto de las arterias de circulacion sobre la accesibilidad podemos hacer un hexbin map para visualizar mejor
montevideo_grid <- st_make_grid(zonas_censales_CODSEC,
n=c(125,125),
what = 'polygons',
square = FALSE,
flat_topped = TRUE) %>%
st_as_sf() %>%
mutate(area = st_area(.)) %>%
mutate(ID = row_number())
ggplot()+
geom_sf(data = montevideo_grid)+
geom_sf(data = zonas_censales_CODSEC,alpha = 0.4)+
theme_bw()
direcciones <- direcciones %>%
st_transform(4326)
direcciones_hex <- st_join(direcciones,montevideo_grid)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
hex_accesibilidad <- direcciones_hex %>%
st_set_geometry(NULL) %>%
group_by(ID) %>%
summarise(Lineas_accesibles = floor(mean(Lineas_accesibles)))
montevideo_grid <- left_join(montevideo_grid,hex_accesibilidad,by = "ID")
montevideo_grid <- montevideo_grid %>%
filter(!is.na(Lineas_accesibles))
ggplot()+
geom_sf(data = zonas_censales_CODSEC,fill = NA)+
geom_sf(data = montevideo_grid,aes(fill=Lineas_accesibles))+
scale_fill_gradientn(colours = c("#7FF2A4","#CAF087","#F6F296","#F2B637","#CD4206","#5F1F03"),
guide = guide_colourbar(barheight = 10))+
labs(title = "Accesibilidad al STM por sector de la Ciudad",
x = "Longitud",
y = "Latitud",
legend = "Accesibilidad")+
theme_bw()
Descargamos el cuadro delimitador mediante la función getbb.
bbox_montevideo <- getbb("Montevideo, Uruguay")
Mostramos el resultado.
bbox_montevideo
## min max
## x -56.43140 -56.02250
## y -34.93806 -34.70185
Descargamos el mapa base mediante la función get_stamenmap.
mapa_montevideo <- get_stamenmap(bbox=bbox_montevideo,
maptype="toner-lite",
zoom=12)
## Source : http://tile.stamen.com/toner-lite/12/1405/2469.png
## Source : http://tile.stamen.com/toner-lite/12/1406/2469.png
## Source : http://tile.stamen.com/toner-lite/12/1407/2469.png
## Source : http://tile.stamen.com/toner-lite/12/1408/2469.png
## Source : http://tile.stamen.com/toner-lite/12/1409/2469.png
## Source : http://tile.stamen.com/toner-lite/12/1410/2469.png
## Source : http://tile.stamen.com/toner-lite/12/1405/2470.png
## Source : http://tile.stamen.com/toner-lite/12/1406/2470.png
## Source : http://tile.stamen.com/toner-lite/12/1407/2470.png
## Source : http://tile.stamen.com/toner-lite/12/1408/2470.png
## Source : http://tile.stamen.com/toner-lite/12/1409/2470.png
## Source : http://tile.stamen.com/toner-lite/12/1410/2470.png
## Source : http://tile.stamen.com/toner-lite/12/1405/2471.png
## Source : http://tile.stamen.com/toner-lite/12/1406/2471.png
## Source : http://tile.stamen.com/toner-lite/12/1407/2471.png
## Source : http://tile.stamen.com/toner-lite/12/1408/2471.png
## Source : http://tile.stamen.com/toner-lite/12/1409/2471.png
## Source : http://tile.stamen.com/toner-lite/12/1410/2471.png
## Source : http://tile.stamen.com/toner-lite/12/1405/2472.png
## Source : http://tile.stamen.com/toner-lite/12/1406/2472.png
## Source : http://tile.stamen.com/toner-lite/12/1407/2472.png
## Source : http://tile.stamen.com/toner-lite/12/1408/2472.png
## Source : http://tile.stamen.com/toner-lite/12/1409/2472.png
## Source : http://tile.stamen.com/toner-lite/12/1410/2472.png
Mostramos el mapa con la función ggmap.
ggmap(mapa_montevideo)
Delimitamos el polígono de la ciudad.
poligono_montevideo <- getbb("Montevideo, Uruguay",
format_out="sf_polygon")
Vemos el resultado.
poligono_montevideo
## Simple feature collection with 2 features and 0 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -56.4314 ymin: -34.93806 xmax: -56.0225 ymax: -34.70185
## Geodetic CRS: WGS 84
## geometry
## 1 POLYGON ((-56.4314 -34.8281...
## 2 POLYGON ((-56.4314 -34.8281...
ggmap(mapa_montevideo)+
geom_sf(data=poligono_montevideo,
fill=NA,
color="darkgreen",
size=1,
inherit.aes=FALSE)+
labs(title="Montevideo",
subtitle="Uruguay",
caption="Fuente: Open Street Map")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Para mantener coherencia con el resto de los trabajos prácticos de la materia, descargamos las calles de la ciudad.
En primer lugar, asignamos el boundig box.
calles_montevideo <- opq(bbox_montevideo) %>%
add_osm_feature(key = "highway",
value = c("motorway",
"trunk",
"primary",
"secondary",
"tertiary",
"unclassified",
"residential"))
Utilizamos la función osmdata_sf para descargar la informació y la mostramos.
calles_montevideo <- osmdata_sf(calles_montevideo)
Nos quedamos con las líneas.
calles_montevideo <- calles_montevideo$osm_lines
Graficamos el resultado.
ggplot()+
geom_sf(data=poligono_montevideo, fill=NA, color="darkgreen", size=1)+
geom_sf(data=calles_montevideo)+
labs(title="Calles",
subtitle="Montevideo, Uruguay",
caption="Fuente: Open Street Map")+
theme_void()
Lo unimos con nuestro mapa base.
ggmap(mapa_montevideo)+
geom_sf(data=poligono_montevideo, fill=NA, color="darkgreen", size=1, inherit.aes = FALSE)+
geom_sf(data=calles_montevideo, inherit.aes = FALSE)+
labs(title="Calles",
subtitle="Montevideo, Uruguay",
caption="Fuente: Open Street Map")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
dim(calles_montevideo)
## [1] 14733 114
calles_montevideo %>%
group_by(highway) %>%
summarise(cantidad=n())
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## Simple feature collection with 6 features and 2 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -56.44871 ymin: -34.93419 xmax: -55.99676 ymax: -34.65824
## Geodetic CRS: WGS 84
## # A tibble: 6 x 3
## highway cantidad geometry
## <chr> <int> <MULTILINESTRING [°]>
## 1 primary 620 ((-56.08889 -34.73053, -56.08974 -34.7317, -56.09039 -3~
## 2 residential 11144 ((-56.03465 -34.84728, -56.03397 -34.84764, -56.03344 -~
## 3 secondary 842 ((-56.11317 -34.72069, -56.11315 -34.72177, -56.11308 -~
## 4 tertiary 1458 ((-56.22758 -34.72739, -56.22695 -34.72737, -56.22586 -~
## 5 trunk 447 ((-56.35509 -34.78715, -56.35661 -34.78021), (-56.10173~
## 6 unclassified 222 ((-56.17819 -34.71774, -56.17826 -34.71765, -56.1795 -3~
Hay 14.733 observaciones, dividiendose en: - 620 primarias. - 11.144 residenciales. - 842 secundarias. - 1.458 terciarias. - 447 troncales. - 222 sin clasificar.
Eliminamos los tramos fuera de nuestra área de influencia.
calles_montevideo <- st_intersection(calles_montevideo, poligono_montevideo)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
dim(calles_montevideo)
## [1] 24056 114
calles_montevideo %>%
group_by(highway) %>%
summarise(cantidad=n())
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## Simple feature collection with 6 features and 2 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -56.43103 ymin: -34.93419 xmax: -56.03084 ymax: -34.70411
## Geodetic CRS: WGS 84
## # A tibble: 6 x 3
## highway cantidad geometry
## <chr> <int> <MULTILINESTRING [°]>
## 1 primary 1126 ((-56.19722 -34.91148, -56.19859 -34.9111), (-56.19859 ~
## 2 residential 17918 ((-56.19148 -34.72162, -56.1917 -34.72148), (-56.18441 ~
## 3 secondary 1640 ((-56.06011 -34.86677, -56.05993 -34.86667), (-56.05993~
## 4 tertiary 2488 ((-56.21423 -34.74944, -56.21422 -34.74923), (-56.21422~
## 5 trunk 692 ((-56.35509 -34.78715, -56.35577 -34.78405), (-56.10173~
## 6 unclassified 192 ((-56.28868 -34.72325, -56.2886 -34.72319), (-56.2886 -~
Hay 14.733 observaciones, dividiendose en: - X primarias. - X residenciales. - X secundarias. - X terciarias. - X troncales. - X sin clasificar.
Eliminamos los tramos fuera de nuestra área de influencia.
Graficamos.
ggmap(mapa_montevideo)+
geom_sf(data=poligono_montevideo, fill=NA, color="darkgreen", size=1, inherit.aes = FALSE)+
geom_sf(data=calles_montevideo, inherit.aes = FALSE)+
labs(title="Calles",
subtitle="Montevideo, Uruguay",
caption="Fuente: Open Street Map")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Descargamos las paradas de omnibus.
En primer lugar, asignamos el boundig box.
paradas_montevideo <- opq(bbox_montevideo) %>%
add_osm_feature(key = "highway",
value = "bus_stop")
Utilizamos la función osmdata_sf para descargar la información
paradas_montevideo <- osmdata_sf(paradas_montevideo)
paradas_montevideo
## Object of class 'osmdata' with:
## $bbox : -34.938056,-56.4313997,-34.7018526,-56.0225006
## $overpass_call : The call submitted to the overpass API
## $meta : metadata including timestamp and version numbers
## $osm_points : 'sf' Simple Features Collection with 204 points
## $osm_lines : NULL
## $osm_polygons : 'sf' Simple Features Collection with 1 polygons
## $osm_multilines : NULL
## $osm_multipolygons : NULL
Nos quedamos con los valores “punto”.
paradas_montevideo <- paradas_montevideo$osm_points
Existen 204 observaciones.
Graficamos.
ggmap(mapa_montevideo)+
geom_sf(data=poligono_montevideo, fill=NA, color="darkgreen", size=1, inherit.aes = FALSE)+
geom_sf(data=paradas_montevideo, inherit.aes = FALSE, aes(color=highway))+
labs(title="Paradas",
subtitle="Montevideo, Uruguay",
color="Tipo",
caption="Fuente: Open Street Map")+
scale_color_manual(values="darkorange")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Se ve que hay puntos por fuera del área de estudio. Se verá primero la dimensión y luego se los quitará con la función st_intersection.
dim(paradas_montevideo)
## [1] 204 35
paradas_montevideo %>%
group_by(highway) %>%
summarise(cantidad=n())
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## Simple feature collection with 2 features and 2 fields
## Geometry type: MULTIPOINT
## Dimension: XY
## Bounding box: xmin: -56.42333 ymin: -34.92366 xmax: -56.02256 ymax: -34.70244
## Geodetic CRS: WGS 84
## # A tibble: 2 x 3
## highway cantidad geometry
## <chr> <int> <MULTIPOINT [°]>
## 1 bus_stop 200 ((-56.42333 -34.75531), (-56.42215 -34.75554), (-56.40782 -~
## 2 <NA> 4 ((-56.22367 -34.72801), (-56.22366 -34.72779), (-56.22364 -~
Vemos que hay 200 paradas y 4 N/A.
paradas_montevideo <- st_intersection(paradas_montevideo, poligono_montevideo)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
ggmap(mapa_montevideo)+
geom_sf(data=poligono_montevideo, fill=NA, color="darkgreen", size=1, inherit.aes = FALSE)+
geom_sf(data=paradas_montevideo, inherit.aes = FALSE, aes(color=highway))+
labs(title="Paradas",
subtitle="Montevideo, Uruguay",
color="Tipo",
caption="Fuente: Open Street Map")+
scale_color_manual(values="darkorange")+
theme_void()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
paradas_montevideo %>%
group_by(highway) %>%
summarise(cantidad=n())
## although coordinates are longitude/latitude, st_union assumes that they are planar
## Simple feature collection with 1 feature and 2 fields
## Geometry type: MULTIPOINT
## Dimension: XY
## Bounding box: xmin: -56.35556 ymin: -34.92366 xmax: -56.03415 ymax: -34.76106
## Geodetic CRS: WGS 84
## # A tibble: 1 x 3
## highway cantidad geometry
## <chr> <int> <MULTIPOINT [°]>
## 1 bus_stop 204 ((-56.35556 -34.82972), (-56.35414 -34.7878), (-56.35386 -3~